%随机改变体系的构型：随机选取一个原子，使其产生一个方向随机、大小为dr *rand的位移
%config0是一个1*N的cell（N为体系的原子数），每个元素是一个1*3数组，记录原子坐标
%limit是空间上的限制，包含三个坐标的上下限，共6个元素
%choice=1时，只移动最后添加的那个原子；choice=0时，除了在原点处的原子外其他所有原子都移动
function [config, en] = change_config(config0,en0,limit ,dr,choice)
    d = 0.88;
    N = length(config0);
    config  = config0;
    en = en0;
    
    ext = 1.5;
    x_min = limit(1) - ext; x_max = limit(2) + ext;
    y_min = limit(3) - ext; y_max = limit(4) + ext;
    z_min = limit(5) - ext; z_max = limit(6) + ext;
    
    if(choice == 1)
        i = N;
        en_i = sim_en_cal(config0, i);
    end
    
    while 1
        if(choice == 0)
            i = ceil(rand*(N-1))+1;
            en_i = sim_en_cal(config0, i);
        end
        list = 1:N;
        list(i) = [];
        r = dr * rand;
        theta = rand * pi;
        phi = rand * (2 * pi);
        dx = r * sin(theta) * cos(phi);
        dy = r * sin(theta) * sin(phi);
        dz = r * cos(theta);
        
        %对新位置进行限制，采用周期性移动的方式
        new = config0{i} + [dx,dy,dz];
        x_new = new(1);y_new = new(2);z_new = new(3);
        if(x_new > x_max)
            new(1) = x_new - (x_max - x_min);
        elseif(x_new < x_min)
            new(1) = x_new + (x_max - x_min);
        end
        if(y_new > y_max)
            new(2) = y_new - (y_max - y_min);
        elseif(y_new < y_min)
            new(2) = y_new + (y_max - y_min);
        end
        if(z_new > z_max)
            new(3) = z_new - (z_max - z_min);
        elseif(z_new < z_min)
            new(3) = z_new + (z_max - z_min);
        end

        judge = 1;

        for j=list
            rj_vector = config0{j} - new;
            rj_square = sum(rj_vector.^2);
            rj = sqrt(rj_square);
            if(rj < d)
                judge = 0;
                break
            end
        end
        
        if(judge == 1)
            config{i} = new;
            en_i_new = sim_en_cal(config, i);
            en = en + en_i_new - en_i;
            break
        else
            
        end
    
    end
    return
end

%计算第order个原子与其他原子的相互作用势之和
function en = sim_en_cal(config, order)
    N = length(config);
    en = 0;
    mlist = 1:N;
    mlist(order) = [];
    for m=mlist
        en = en + LJ_pot(config{m},config{order});
    end
    return
end

%取原子直径为0.88sigma，截断半径为3.0sigma
%pt1,pt2均为1*3的数组，记录原子中心在三维直角坐标系下的3个坐标分量
function pot = LJ_pot(pt1, pt2)
    %根据pt1,pt2两点的坐标计算两点距离
    rvector = pt1 - pt2;
    r_square = sum(rvector.^2);
    r = sqrt(r_square);
    pot = r^(-12)-r^(-6); 
    
%     if(r < 0.88)
%         pot = 0.88^(-12)-0.88^(-6);
%     elseif(r > 3)
%         pot = 0;
%     else
%         pot = r^(-12)-r^(-6);   
%     end

    return
end